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Abstract 

A bilinear quadrature numerically evaluates a continuous bilinear map, such as 
the L 2 inner product, on continuous / and g belonging to known finite-dimensional 
function spaces. Such maps arise in Galerkin methods for differential and integral 
equations. The construction of bilinear quadratures over arbitrary domains in R d is 
presented. In one dimension, integration rules of this type include Gaussian quadrature 
for polynomials and the trapezoidal rule for trigonometric polynomials as special cases. 
A numerical procedure for constructing bilinear quadratures is developed and validated. 


1 Introduction 

Classical quadratures such as Gaussian and trapezoidal rules accurately evaluate continuous 
linear functionals such as 

/ f(x)w(x) dx 
Jn 

for / in a finite-dimensional space of continuous functions. Bilinear quadratures evaluate 
continuous bilinear forms such as the weighted L 2 inner product 


(f,g)L 2 = / f{x)g(x)w{x)dx 
Jn 


or the weighted H 1 inner product 


(/>s)m = / X 


( df . . dg \ 
I- aij(x) -I 


' i,j =1 


\dx 


dxjJ 


+ f{x)g(x) dx 


on finite-dimensional spaces of continuous functions /, g on 17 C R d . 

L 2 inner products compute orthogonal projections onto subspaces, while H 1 inner prod¬ 
ucts provide local solutions to elliptic problems, a key ingredient of the finite element 
method. For example, let 9 C R d be a smooth bounded domain, let L be a uniformly 
elliptic operator, / £ L 2 (Q), g £ L 2 (9n), and 7 £ L°°{d f2). Consider the Robin problem 


Find u £ satisfying 

Lu = f in Cl, ju + = g on 


( 1 . 1 ) 


When L = —A, then if a bilinear form a : H 1 ( f!) x —>• R is defined by a(u,v) = 

f n Du ■ Dv + fgQ'Juv, the weak formulation to (1.1) seeks u £ i7 1 (fl) satisfying 


a(u,v) = {f,v) L 2 (n) + (g,v) L 2 for all v £ H ^fi). 


( 1 . 2 ) 


The Galerkin method constructs an approximate solution to (1.2) by choosing finite-dimensional 
function spaces Do, Go and seeking uq £ Do satisfying 


a(u 0 ,u 0 ) = (f,v 0 )mn) + (g,v o)L2(3fi) for all v 0 G Go- 


(1.3) 
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The linear system (1.3) is solved in a basis, which requires computing a number of L 2 inner 
product integrals. These integrals should be computed both efficiently and accurately. 

Efficiency is achieved by using the fewest function evaluations possible. When d > 1, 
the optimal efficiency of a classical quadrature is unknown. For a bilinear quadrature, the 
minimum number of function evaluations is equal to the dimension of function space being 
integrated. The inner product of two functions /, g belonging to given finite-dimensional 
function spaces is computed by the formula 

(f,g) = f(x)*Wg( y), (1.4) 

where /(x) £ R m and g(y) £ R” are evaluations of / and g at sets of points x and y in 
Q, respectively, and W is a matrix. The rank of the bilinear form is equal to the rank of 
W, hence the minimal number of required function evaluations is equal to the dimension of 
that function space. 

Accuracy is achieved by defining and minimizing integration error. In a bilinear quadra¬ 
ture, this is a nonlinear optimization problem for x, y, and W in (1.4), and is solved using a 
Newton method for an appropriate objective function [CRY99, BGR10, XG10]. In this pa¬ 
per an objective function is developed and demonstrated to yield numerically useful bilinear 
quadrature rules in a general setting. 

Numerical evaluation of inner product integrals has been studied in [BD71, McG79, 
Gri80, BGR10, Chel2] and as “bilinear quadrature” in [LZ87, Kno07]. This paper bor¬ 
rows some of the framework from these past works but develops and utilizes a different 
optimization procedure to produce quadrature rules. 


2 Theory 

2.1 Abstract formulation 

In this section the problem of evaluating a general continuous bilinear form on a pair of 
Banach spaces is considered. Results are given in great generality so that they apply to any 
continuous bilinear forms. Later, these results are applied to useful special cases such as the 
L 2 and H 1 inner products. 

Definition 2.1. Let T and Q be real Banach spaces. Then a bilinear quadrature of order 
(m, n) on T x Q is a bilinear form Q defined by linear maps L\ : T —> R m and L 2 : G —t R n 
and a bilinear map B : R m x R n — > R, such that, for each / £ T and g £ G, 

Q{f,g) = B(L 1 f,L 2 g). 

Definition 2.2. Let A, G be real Banach spaces with a continuous bilinear form (•,•) : 
T x G —> R. Finite-dimensional subspaces JoC J and Go C G are a dual pair if 

V/€J„\ {0}, 3g £ Go such that (/, g) ± 0, 

Vs G Go \ {0}, 3/ e To such that (/, g) ± 0. 

If To, Go are a dual pair then dim(Jx)) = dim(t/o)- 

Definition 2.3. Let T, G be real Banach spaces with a continuous bilinear form (•,•) : 
T x G —t R, and let To d *F and Q§ d Q be a dual pair. A bilinear Quadrature Q on x Q 
is exact with respect to To x Qo if 

(f,9) = Q{f,9) for every f £T 0 ,g £ Go- 
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Such a bilinear quadrature evaluates the bilinear form on To x Go exactly. We have the 
diagram 

X t/’o -»•—H ; R m X R n 



B(-r) 

> ' 

R 


If the parent spaces J r , G are implied, we will abuse notation by referring to an exact bilinear 
quadrature on To x Go- 

Remark. If T , G are infinite-dimensional and Q is exact on J- 0 x Go, then 

sup \Q(f,g) - {f,g)\ = oo, 

/e^.seC 

so a bilinear quadrature can only be accurate on finite-dimensional subspaces. 

Lemma 2.4. Let Jo C J and Go C G be a dual pair, and let L\ : T —» R m , L 2 : £ —> R n be 
linear. Then there exists bilinear B : M m xl" — > R such that the map Q(f , g ) = B(L\f, L 2 g) 
is an exact quadrature on To x C/o */ and only if Li|_f 0 and L- 2 \g o are both injective. 

Proof. Suppose B exists. If /, / € To are distinct then there exists g £ Go such that 

0 ?(f- f, 9) = Q{f - f, 9) = B(Li (/ - f),L 2 g), 

so L\f Lif and Li|jr 0 is injective. Similarly for L 2 \g o . 

Suppose Li\jr a and L 2 \g Q are injective. Their Moore-Penrose pseudoinverses (Li| jro ) + 
and (L 2 \g o ) + left-invert Li and L 2 , respectively. Define a bilinear map on R m x R” by 

B(x,y) = ((ii|jr 0 ) + x, (L 2 \g o ) + y'j . 

Then, for all / G To, g € Go, 

B{L\f, L 2 g) = <(Li|jr 0 ) + Ji|jF 0 /, (i 2 | go ) + J 2 | eo <?) 

= (/>5>- 


□ 

From Lemma 2.4 a necessary condition for an exact bilinear quadrature is that to > 
dim To , n > dimt/o. Minimal order is achieved when m = dim To, n = dimC/o and B{x,y) 
is uniquely given by 

B(x,y) = ( s {L\\ :Fa )~ 1 x, (L 2 \g o )~ l y) . 

Exact bilinear quadratures are not unique, as there are many possible linear maps L\,L 2 . 
Furthermore, B may not be unique, since if n > dim(J' 0 ), then Li\jr a has infinitely many 
left inverses. Therefore, a method is needed to choose among the infinitely many bilinear 
quadratures. One metric of quality is that, in addition to its exactness on To x Go, the 
bilinear quadrature also approximates (f,g) for some set of g' s outside of Go- 

Definition 2.5. Let To C T and Go C G be a dual pair, and let Gi C G be another 
finite-dimensional subspace such that 

Gi C To := {g G G : ( f,g ) = 0 for all / € J'o}- 

Let Q be a set of bilinear quadratures exact on To x Go- Then Q € Q is called minimal on 

Gi if 

Q = argmin a(Q;T 0 , Gi), (2.1) 

QeC 
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where 


a(Q]T,o,Gi)'■= max 


\Q(f,g)\ 


o^feJ^o 

If Q is minimal on Gi, then it approximates the pairing of J~q and Go © Gi- Precisely, if 
/ € Jo, 9 € Go © Gi, and we write g = g 0 + gi with < 7 , £ Gi, then 


I Q(f,g) - (/, g)\ = \Q(f,gi)\ < <x(Q;F 0 ,Gi)\\f\\M\ g . 


( 2 . 2 ) 


Thus, minimizing a(Q; Fo,Gi) will improve the approximation. 

One important special case for bilinear quadratures is the symmetric case, which is when 
F = G is an inner product space. In this case, a bilinear quadrature computes an orthogonal 
projection. 

Definition 2.6. Let Jo and Go be a dual pair in an inner product space. Let {/)} be an 
orthonormal basis for Jo. Given a bilinear quadrature Q exact on Jo x Go, the approximate 
orthogonal projection onto Jo arising from Q is the linear map Pq given by 

p q{9 ) = ^2Q{fi,g)U 


An error estimate for orthogonal projections similar to (2.2) is given later in Theorem 2.7. 


2.2 Integral formulation 

In this section, the bilinear quadrature framework is applied to the evaluation of Sobolev 
inner products on function spaces. Let O C R d be a bounded domain, and let J = G = 
C r (il), r a non-negative integer, equipped with a Sobolev inner product 


(f,9)H- = E (D a f,D a g) L 2 {n) 

\at\<s 

for s < r. 

Choose a dual pair Jo, Go in J- Exactness on Jo x Go requires linear maps L\ : C r (fl) —> 
U, L 2 : C r (n) —► U, and bilinear form B : R m x R n —► R so that for every / £ Jo, g £ Go, 
B(L 1 f,L 2 g) = ( f,g)H° ■ 

Appropriate linear maps Li, L 2 are pointwise evaluations at particular points in 12. Thus, 
for the points x = (aq,..., x m ) £ f2 m , y = (t/i,..., y n ) £ fl n , define 


70i)" 


3(yi)" 


, L 2 g := g( y) = 




_gi.Vn)_ 


Given bases /3 = {/ 1 ,..., /*,} for J 0 and {< 71 ,..., g^} for Go, let M £ M fcxfc be the Gram 
matrix with entries 

Mij {fi, 9j) H s • 

Since Jo, Go are a dual pair, M is invertible. Define matrix functions 



Fi(xi) . 

• fk(x 1)" 


’Gi(yi) ■ 

■ 9k{yi) 

J(x):= 



, G( y):= 




1 (*^m) 



_Gi(y n ) ■ 

■ 9 k(y n )_ 


To make L 1 and L 2 are injective, choose x,y, such that J(x) and G( y) have full column 
rank. If B{y,w) =v*Ww for all v,w for an m x n matrix W, then the bilinear quadrature 
is exact if and only if 

J(x)*ITG(y) = M. (2.3) 
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Therefore a bilinear quadrature rule 


Q(f,g) = f(x)*Wg{ y) (2.4) 

evaluates {f,g)H s exactly for any / C Jo ,g £ Qq. The corresponding approximate orthogo¬ 
nal projection onto Jo is 

k 

P Q (g) = J2lMx-T w 9( y)]U 

In the basis /3, the approximate projection is computed by 

[P Q { 9 )b=F^)*Wg( y)eM fc . (2.5) 


Good values for the matrix W and evaluation points x, y must be determined. Without 
loss of generality, suppose that the bases {/*} and {gj} are iJ-orthonormal in G r (fl). Select 
finite-dimensional Q\ C G r (fl) for the minimization (2.1) and define the feasible set Q to 
be all quadratures of the form (2.4) satisfying (2.3). If { 71 ,... , y p } is an orthonormal basis 
for Qi, define 

’71(^1) ••• 7 P (xiY 


T(x):= 


<E K” xp . 


ji(x n ) ... 7 p {x n ) 


Then (2.1) can be reformulated as 


min cr(Q; J 0 ,£i) = min max \Q{f,g)\ 

Q&Q QeC 3661 ,||s||g=i 

/e J),||/||f=i 

= min max | 6 *F(x)*WT(y)a| 
x <yW a ,beR fc 

||«|| 2 =|| 6 || 2=1 

= min fj 1 (F(x)*WT(y)) subject to f(x)*lbG(y) = M, (2.6) 

x,y .W 


where cri(A) is the leading singular value of a matrix A. Minimization (2.6) is independent 
of x, since by (2.3) F(x.)*W = ML, where L is a left inverse of G(y). Therefore x is chosen 
by performing a similar minimization on the left, setting an orthonormal basis {} for a 
space T\ C So - , defining the corresponding matrix function A(x), and minimizing 

min ( 71 (A(x)*WG(y)) subject to F(x)*WG( y) = M, (2.7) 

x,y,W 

where similarly the dependence of (2.7) on y may be dropped since WG{ y) is equal to L* M, 
where L is a left inverse of F(x). 

In the symmetric case Jo = So, M = J, F\ = Si, and m = n = k with x = y, 
minimizations (2.6) and (2.7) are equivalent and simplify to 

mincri(J(x) _ 1 r(x)). ( 2 . 8 ) 

X 


In subsequent sections special attention is given to the symmetric case because it is used 
for evaluating orthogonal projections. 


2.3 Error estimates 

In this section, upper bounds on several error quantities in computing an approximate 
orthogonal projection of the form (2.5) are estimated. 

Theorem 2.7 (Euclidean norm error estimate). Let Jo, So be a dual pair in an inner 
product space J and Q a bilinear quadrature of the form (2.4) that is exact on Jo x Qq. 
Let Pq be the approximate orthogonal projection onto Jo arising from Q with coordinate 
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representation (2.5). If P is the exact orthogonal projection operator onto To, Q\ C Jyj 1 , 
and g = g 0 + gi £ Go © Qi such that gt £ Gi, 


\\[P Q (g)-P(gM 2 <cr(Q;To,Gi)\\gil 

where || - ||2 is the Euclidean norm. 

Proof. This is essentially the same as (2.1). Since (fi,g) = Q(fi,go ), then 


1/2 


\\[Pq(9) - P(g)b II 2 = EI«(/<.») - </»s>f 


\i=l 
' k 


1/2 


= ElQ(/i,s)-Q(/i,sb)l s 

\i=l 

/ k \ 1/2 

= (ElWorfJ 


1 ' 

= max— r r '52a i Q{f i ,g 1 ), 

a ^° IMh 

where a = (af) £ M fc . Each f £ To can be written as / = JU a ifi > so 

k 


max t:— r- CLiQ(fi,gi)= max 

a 2 ^ OA/e^-o 

2=1 


Q(/,Si) 


<ff(Q;^b,ai)||ffi||. 


(2.9) 


0 


Theorem 2.7 provides an error bound for an approximate orthogonal projection when 
the projected function g is in Go®Gi- If To is a space of polynomials, then it is also useful 
to obtain an error estimate that depends on the regularity of g. 

Theorem 2.8 (Uniform norm error estimates for polynomials). Let T = 0(12) with U C 
a hounded , convex domain, equipped with the L 2 inner product. Let Tq be the set of 
multivariate polynomials of degree at most n with an orthonormal basis ft = {/$}, let P : 
T —► T be the orthogonal projection onto To, and suppose Pq is an approximate orthogonal 
projection onto Tq with coordinate representation (2.5). There exist a constant C > 0 such 
that for every g £ C n+l {Vt), then 

II [Pg - P Q gb Hoc < c\\D n+1 g\\ LBB) (2.10) 

where 

P" +1 <?lk» := E max \D a g(x) |. 

|a|=n+l 

Proof. Using (2.5) and the exactness of Pq on Tq, then writing g = Pg+ (/ — P)g = go + gi, 
we have 


I \[Pg - PQg\p\\oo = \\F*Wgi (xJHoo 

< ||0*IU|| 00 —>-oo \\(i-P)g\\ c° 

< ||E*IU||oo^oo||(/ - P){g - q)Wc°, 

where q is any element of Tq and || • ||oo->-oo is the induced matrix 00 norm. Then 
l|[P<7 - P Q g\h lloo < ||E*IU||oo^oo(l + ||P||co^ C o)||3 - q\\co, 
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where the C° operator norm of P is given by 


\\P II 


c°- 


>c o = max 

x€:Q, 




dt. 


By the Deny-Lions/Bramble-Hilbert lemma [EG04], for all g £ C' rl+1 (fl) there exists a 
constant Cbh > 0 (dependent on n and fl) such that 


inf 

q&Fa 


\\g - q\\c° < C B H\\D n+ g\\ L 


which combined with the previous inequality yields the desired result with 


C = ||F*W|| 00 ->. 0 o(l + \\P\\c°^c°)Cbh- 


□ 

In the presence of round-off error in function evaluation, the conditioning of an approx¬ 
imate orthogonal projection is also important to quantify. 

Theorem 2.9. Let Pq be an approximate orthogonal projection of the form (2.5). IfSgfy) 
is the absolute error in computing g(y) and SPQ(g) is the resulting projection absolute error, 
then with respect to a vector norm || • ||, 

1MM< ll*g(y)ll 

\\[Pq(9M - ll<?(y)ll ’ 

where k = n(F* (x)W) is the matrix condition number with respect to || • ||. 

2.4 Classical and bilinear quadratures on univariate polynomials 

In this section we review Gaussian quadratures and show they are a special case of a bilinear 
quadrature in one dimension. We then propose a way to generalize to quadratures evaluating 
inner products of polynomials on multidimensional domains. 

Definition 2.10. Let Cl C be a connected domain. A classical quadrature q of order n 
on Cl is a linear functional defined by a set x = (x\, ... , x n ), Xi £ Cl, called the nodes , and 
a vector w £ M", whose components are called the weights, such that for any / £ C(Cl), 

n 

<?(/) = w*/(x) = w if( x i)- 
1=1 

Furthermore, if Tq is a subspace of C(Cl) and p, is a Borel measure, then q is said to be 
exact on Pq if 

?(/) =[ f dp 
Jn 

for all f £ Pq. 

Let P n be the space of univariate polynomials of degree up to n , I an open interval, and 
p a finite absolutely continuous Borel measure on /. 

Definition 2.11. Suppose P 2 ra -i is /i-integrable on I. Then a Gaussian quadrature of order 
n on I is a classical quadrature of order n on I that is exact on P 2 n-i with respect to p. 

The advantages and disadvantages of the theory of quadratures for polynomials are 
rooted in existence and uniqueness result for Gaussian quadratures. 

Theorem 2.12. Suppose P 2 n-i is p-integrable on I, and let {4>k} denote any set of L 2 (I, p)- 
orthonormal polynomials such that deg(</>&) = k. Then the following sets are equal: 
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1. The zeros of <f> n . 

2. The eigenvalues of the symmetric bilinear form on P n _i given by 

B{f,g) '■= J xf (x)g(x) dfi = {xf(x),g(x)) L 2 ^ I ^ ) ■ 

3. The nodes {xi} of a Gaussian quadrature of order n on I. 

Proof. (1 •<==> 2) Since B(-,-) is symmetric it is diagonalizable with n real eigenvalues {A.;}. 
If a polynomial ipi{x) is an eigenvector for A,, then for 0 < j < n — 1, 

(■.xil)i{x),<t>j{x )} = A i(ipi(x),<i>j(x)} ==> {(x - \i)ifi(x),(f)j{x)) =0. 

Since (x — A i)ipi(x) € P n for each i, and the only polynomials in P„ that are orthogonal to 
each of <j>o, ■■■, 4>n- 1 are multiples of <j> n , then each (x — A*) is a factor of <j> n (x). Thus <j> n (x) 
is a multiple of (x — Ai )... (x — A n ) and its zeros are the eigenvalues of £?(•, ■). 

(2 <*==>- 3) Suppose a Gaussian quadrature with weights and nodes {a^} exists. 
With respect to the basis of orthonormal polynomials {4> k }, the bilinear form B(-,-) has a 
symmetric matrix represention B with entries given by 

/ n 

x(fi(x)<j)j(x) dg, = Yw k x k Mxk)Mx k ). 

fe=l 


If 



y/WkMxk) 


Xi 

Uk = 


,x = 



y/Wktyn— 1 (Xk) _ 


x n_ 


then 

n 

B = Y x kUkU* k = UXU*. (2.11) 

fc=i 

Since Sij = w k<t>i(xk)4 > j{x k ), then / = UU* and U is a unitary matrix. Then (2.11) 

is the unitary diagonalization of the symmetric matrix B with eigenvalues given by the 
x k ’s. □ 

Remarks. Theorem 2.12 shows that if a Gaussian quadrature of order n exists, its nodes are 
the zeros of f> n . The existence proof is completed by showing the weights exist and satisfy 

n 

1 M = ■ 

3 =0 

Therefore, taking the square root ^ydik is legitimate [DR84]. Theorem 2.12 also provides an 
efficient method to construct these quadratures. The matrix B in (2.11) is tridiagonal, so 
its eigenvalues can be calculated quickly, even for very large n [GW69]. 

Gaussian quadrature is optimal for integrating polynomials on an interval, but does not 
extend readily to higher-dimensional domains. The zeros of a multivariate polynomial are 
generally not isolated (consider f(x,y) = xy) so they cannot all be used as nodes of a 
classical quadrature. Additionally, the connection between nodes and eigenvalues no longer 
holds since the eigenvalues are only scalars (the connection extends to two-dimensional 
domains with complex eigenvalues [VR14]). 

Bilinear quadratures make sense in any dimension, yet contain Gaussian quadrature as 
a special case. Consider a classical quadrature with nodes x = {.t,} and weights w = {w,}. 
If a function h(x) = f(x)g(x) with /, g belonging to function spaces J 7 , Q respectively, then 
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the classical quadrature q evaluated on h is the same as a bilinear quadrature Q on /, g 
given by 

W\ 


Q{f,g) = f{*Y 


g(x) = w *h{x) = q(h). 


w 


n 


The matrix W is diagonal with entries given by the weights of the classical quadrature. 
Thus for a general bilinear quadrature of the form (2.4) the entries of W can be viewed as 
analogues of the weights. 


Theorem 2.13. The nodes of a Gaussian quadrature of order n are the same as the points 
x in the unique bilinear quadrature of order (n,n) on P n _i x P n -i that is minimal on 
span{0 n }, where <f> n is the orthonormal polynomial of degree n. Furthermore, the matrix 
W in (2.4) is a diagonal matrix whose diagonal entries are the weights of the Gaussian 
quadrature. 


Proof Let 0o,..., (f> n -i be the orthonormal polynomials up to degree n — 1 such that 
degcfk = k, x = (xi,..., x n ) G 11 ", and define 


$(x) = 


</>o(zi) ■ ■ 


00 (Xn) ■ ■ 

Then the minimization problem (2.8) becomes 


min (j i dhx) 
xGQ™ 1 


4>n- l(*l) 
0n—l(^n) 

4>n(x l) 

0n (^n) 



This is uniquely minimized (up to reordering of the xf s) when x is the set of zeros of 0 ra 
in which case a i = 0. The corresponding bilinear quadrature Q exactly evaluates products 
where one polynomial has degree n — 1 and the other has degree n. Then Q has the same 
evaluation points as a bilinear quadrature formed from the Gaussian quadrature. Since W 
is unique, then it must be equal to the diagonal matrix with entries given by the weights of 
the Gaussian quadrature. □ 

The above result suggests that a good way to accurately compute inner products of poly¬ 
nomials on a multidimensional domain is to utilize a symmetric bilinear quadrature that is 
exact on P„ x P„ and minimal on P„+i ClP^. Just as Gaussian quadratures accurately inte¬ 
grate nearly polynomial functions accurately, bilinear quadratures constructed in the above 
manner are expected to evaluate inner products of nearly polynomial functions accurately. 
Numerical results for these quadrature are shown in section 3. 


2.5 Classical and bilinear quadratures on trigonometric polynomi¬ 
als 

For the space of trigonometric polynomials 

T n _i = spanjl, sinx,..., sin (n — l)x, cosx,..., cos (n — l)x}, 


it is known that the (n + l)-point trapezoidal rule 


Tra(p) := -p(0) + -p{2n) + — ^ p f 
n n n ' V 

.7=1 x 


2ttj 

n 
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is exact for integrating all p £ T„_i over the interval [0, 27 t]. Since T„_i is a rotationally- 
invariant function space on the circle R/27tZ, the trapezoidal rule yields a family of n-point 
classical quadratures for T„_i given by 


r 2 ir If_^ 27f 

/ p(x) dx = — ’S~~]p(xj) for all p £ T n _i, Xj+ 1 — £ j = —. 

Jo n ]=0 n 


( 2 . 12 ) 


When n is odd, the above trapezoidal rule quadrature is a special case of a bilinear quadra¬ 
ture: 


Theorem 2.14. Let n > 0 be an odd integer. Then the set of classical quadratures on T n _i 
in (2.12) are equivalent to symmetric bilinear quadratures of order ( n,n ) on Tr n ~i )/2 x 
T( n -i)/ 2 that are minimal on spanjsinnx, cosnx}. 


Proof. Set n = 2k + 1. Since T n _i is rotationally invariant on the circle, then if Q is 
a symmetric bilinear quadrature on Tj~ x of the form (2.4), a(Q) is invariant under 
rotations of the evaluation points x = (aq,..., x n ). Therefore without loss of generality 
X\ = 0 ,Xj £ [ 0 , 27 t ). Define 


F(x) 



1 

1 


1 

1 

e ~ikx 2 

p ikx 2 


e ~i(k+l)x 2 

e i(k+l)x 2 

g—ikx n 

gikx n 

e ~i(k+l)x n 

e i(k+l)x n 


and it suffices to prove that choosing Xj = 27 t j/n solves the minimization problem ( 2 . 1 ). 

If xj = 2/Tj/n, then the first column of F(x) is the second column of T(x), and the last 
column of F(x) is the first column of T(x). Therefore, in this case, J^(x) _ 1 r(x) = \e n ei], 
where e ? - is the j-th standard coordinate vector, hence cri(F(x) _ 1 r(x)) = 1 . 

We claim that for any choice of nodes x £ [0, 27r)" with X\ = 0, 

a 1 (F(x)" 1 r(x)) > 1. (2.13) 


If (2.13) is established, then setting x 7 = 2n(j — 1 )/n yields a minimal quadrature in Q. 
To show this, let u be the first column of F(x) _ 1 r(x). We will show that ||u ||2 > 1, from 
which (2.13) follows. The column u satisfies the equation 


F(x)u 


1 


1 

a —i(k+ l)x 2 


0 —i(k+ l)x„ 


which is equivalent to the Vandermonde system 


1 1 

1 


1 

1 e “ 2 . 

e z(n- l)x 2 

u = 

e~ix 2 

1 e“" . 

e i(n—l)xn 




Setting Zj = e lXj , then the entries of u are the coefficients of a degree n — 1 complex 
polynomial p(z) such that p(zj) = 1/zj. Setting q(z) = zp{z ), it suffices to find a degree n 
polynomial q(z) such that g(0) = 0 and q(zj) = 1. Such a q is unique and 

n 

q(z) = 1 - II _ z / z i) ■ 

3 =1 

Then the leading coefficient of q , which is also the leading coefficient of p, has absolute value 
1 , and hence ||u|| 2 > 1 . □ 
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Remark. The trapezoidal rule uniquely generates a minimal bilinear quadrature, since in 
that case q(z) = az n for some |a| = 1. If Xj 's are not equispaced, q(z ) has some nonzero 
lower-order coefficients. 

2.6 Lobatto quadrature and the non-invertible case 

While minimizing the number of evaluation points will reduce the cost of evaluating a 
quadrature, it may be advantageous to use more points than is optimal in order to improve 
accuracy. One example is Lobatto quadratures for polynomials of one variable, which use 
more points than Gaussian quadratures. In this section, we observe Lobatto quadratures 
are a special case of a bilinear quadrature where extra evaluation points are used, in which 
case the matrix function F(x) is non-invertible. The formulation of Lobatto-like bilinear 
quadratures on general domains is given. 

Definition 2.15. Let P 2n _i be ^-integrable on an interval / = [a, b\. Then the correspond¬ 
ing Lobatto quadrature is a classical quadrature of order n + 1 exact on P 2n _i with respect 
to g such that if xq, ■ ■ ■ , x n are the nodes, then Xq = a and x n = b. 

Theorem 2.16. Suppose P 2n -i is p-integrable on I, and let {</>&} denote the unique set 
of orthonormal polynomials such that deg (<j>k) = k. Then there exists a unique Lobatto 
quadrature of order n + 1, and the interior nodes {xi : 1 < i < n — 1} are the zeros of 

A Lobatto quadrature of order n + 1 corresponds to a symmetric bilinear quadrature 
that is exact on P„_i x P n _i and minimal on P ra in which the W matrix is diagonal and 
the matrix F(x) is given by 


<t>i{a) 

M a ) 


M x i) 

4>i{x n - 1 ) 

(Xn—l 

Mb) 

4>k(b ) 


Unlike in the Gaussian quadrature case, the matrix F = F(x) is not square, so there 
exists infinitely many matrices W satisfying (2.3). Therefore, the simplified minimization 
condition ( 2 . 8 ) cannot be employed, and one must optimize over both the quadrature nodes 
x and matrices W. In general, suppose F is m x k and G is n x k : both with full column 
rank. Then all matrices W satisfying (2.3) are of the form 

W = (F*)+MG + + Y - FF+YGG+, (2.14) 

where Y is an arbitrary mxn matrix. In the symmetric case F = G and M = /, a minimal 
bilinear quadrature is found through the unconstrained minimization 

Find Y e ]R mxm and x minimizing a x ( F + T + F*Y(I - FF + ) V) (2.15) 

While computationally more expensive, this optimization procedure can be used to compute 
symmetric Lobatto-like bilinear quadratures. First fix points Xq that the bilinear quadrature 
is required to use, then construct the (typically non-square) matrix function F(x o, x), where 
only the points x are varying. Then minimize according to (2.15). This procedure is 
applicable for arbitrary domains fl, any space of continuous functions, and any inner product 
on that space. 

2.7 Change of variables 

For a bilinear quadrature computing an L 2 (fl) inner products, a bilinear quadrature can be 
cheaply constructed for L 2 (<f>(fl)) inner products, where $ is an affine invertible change of 


11 of 20 





Christopher A. Wong 


variables. For continuous functions /) on fi, set 

= fi( x )\ det(D < f>)| _1/2 . 

Then 

{fi,fj)L^(^(n))= / fi{y)fj(y)dy = / fi(x)fj{x)dx= 

2) J n 

The Jacobian £)$ is constant when $ is affine, so if the bilinear quadrature on L 2 (fl) exact 
on Jo x g 0 is 

Q(f,g) = f(x)*Wg(y), 

a bilinear quadrature on L 2 ($(fl)) for (J 0 o (h” 1 ) x ( Q 0 o $ _1 ) is given by 

Q(f, 9 ) = fmx))*Wg($( y)), W = W\ det( D*)\-\ (2.16) 

For an H s inner product with s > 0, in general a new bilinear quadrature cannot be 
cheaply constructed under a change of variables. However, when <h(x) = A Ux + b is affine 
with A G R and U a unitary matrix, a change of variables can still be performed at low 
cost. Let W be the matrix in a bilinear quadrature of form (2.4) computing H 1 ( f2) inner 
products. Then write W = Wq + W±, where Wo is th e matrix for a bilinear quadrature that 
computes L 2 (Cl ) inner products. Then a new bilinear quadrature for H ^^(fl)) is formed 
with matrix 

W= |A|“ 1 Wo + |A|- 3 Wi 

and evaluation points mapped by <F 

3 Computation 

In this section, a basic numerical procedure to produce symmetric bilinear quadrature rules 
is described. Afterward, some numerical examples of bilinear quadrature rules are presented. 

3.1 Orthogonalization 

For a function space Jo, one may initially have a numerical routine to evaluate (up to 
machine precision) basis functions ipi, ... for Jo that are not orthonormal. Assuming 
that the inner products = Mij can be computed exactly, J(x) is computed from 

T(x) and Gram matrix M by 

1. Compute the lower triangular matrix L in the Cholesky factorization M = LL*. 

2. For a given x, perform a lower-triangular solve on the matrix equation *F(x)* = LZ. 

3. Set J(x) = Z*. 

The same procedure can be used to produce an orthonormal basis for the function space Ji 
that the bilinear quadrature is minimized against. 

3.2 Nonlinear optimization 

For the invertible symmetric case we have reduced our problem to the minimization problem 

( 2 . 6 ): 


Find x minimizing a\ (J(x) 1 T(x)) . 

This is a nonlinear optimization problem in d ■ k variables, where d is the dimension of the 
integration region Q and k = dim(Jo). 

The problem of minimizing the largest singular value of a matrix function A(x) is equiv¬ 
alent to minimizing the largest eigenvalue of the symmetric positive semidefinite matrix 


12 of 20 



Christopher A. Wong 


A*(x)A(x). This type of the eigenvalue optimization problem has been extensively studied 
in its own right; see [OW95] [SF95]. 

Often -F(x) and T(x), but not their derivatives, can be accurately computed. Also, the 
multiplicity of the largest singular values are generally unknown. Consequently, a quasi- 
Newton method is ideal for the optimization procedure. The objective function is non- 
convex and typically has multiple local minima, so the optimization procedure is run with 
many initial guesses. Furthermore, in the presence of many nearby local minima, after 
each convergent result, the computed points can be perturbed by a small value <5 and the 
procedure run again with perturbed points as another initial guess. This is repeated until 
suitable convergence. While this procedure may be expensive, computing a quadrature 
is typically a one-time cost, after which the quadrature can be used repeatedly for its 
applications. 

In our numerical experiments, we employ a quasi-Newton method with BFGS updates 
as implemented as part of Matlab’s fminunc routine [Bro70, Fle70, Gol70, Sha70]. Since 
F(x)- 1 r(x) is a small, dense matrix, its norm is computed by calculating its full SVD. 
Up to 10 5 initial random points uniformly distributed across the domain are used, and the 
procedure is iterated until convergence in double-precision arithmetic. 

For our numerical implementation we do not reinforce the constraint that the evaluation 
points Xi remain in the integration domain Q. While in general the full constrained mini¬ 
mization problem may be necessary, we have empirically observed that it is not necessary 
for quadratures on polynomials. This can be explained by observing that the orthogonal 
polynomials grow rapidly outside of U; thus points outside the domain are not expected to 
be good candidates for the solution to the minimization problem. 

Remarks. In the case of polynomials it is possible to accurately compute the gradients of 
U(x) and F(x), in which case a quasi-Newton method may be unnecessary. The BFGS 
method has been chosen since it is robust for different function spaces. 

3.3 Bilinear quadratures on triangular domains 

In practical applications one of the most important cases to consider is the L 2 product 
of polynomials on a simplex. For example, in the finite element method one typically 
solves a two-dimensional PDE locally on polynomials supported on triangular domains. The 
discretization requires computing a number of inner products. In this section we compute 
bilinear quadratures that are exact on polynomials on a triangular domain. 

Because the space of polynomials is affine-invariant it suffices to find evaluation points 
for polynomials on a reference triangle. Given a bilinear quadrature on a reference triangle a 
bilinear quadrature for polynomials on any other triangle can be cheaply obtained using the 
change of variables formula (2.16). A basis of orthogonal polynomials on the right triangle 
with vertices (—1, —1), (—1,1), (1, —1) is given by 

KmA^V) = (^) m Pm ( 2g: +^ + 1 ) P^+l,0 {y)j (3.1.) 

where P m is the mth Legendre polynomial and P“ ,/3 is the nth Jacobi polynomial with 
parameters a,/3. These functions can be computed efficiently and stably as in [XG10]. 

Using this basis, symmetric bilinear quadratures exact for the L 2 inner product over this 
right triangle on J-q = P„ and minimal on T\ = P n +i n P^ were computed. The minimal 
number of evaluation points were used, in which case the number of points required is 

k = dim(P„) = ^ . 

In Table 1, for each computed bilinear quadrature rule, the minimized largest singular value 
a = cri(U(x)- 1 r(x)) is given, along with the oo-norm condition number of the matrix for 
the approximate orthogonal projection. 
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n 

k 

a 

Koo (F*W) 

0 

1 

0.00000 

1 .00000e+0 

1 

3 

0.14507 

2.82218e+0 

2 

6 

0.30373 

6.29185e+0 

3 

10 

0.47762 

1.15455e+l 

4 

15 

0.65817 

2.03810e+l 

5 

21 

0.78394 

3.39955e+l 

6 

28 

0.87930 

4.71065e+l 

7 

36 

0.95305 

8.48889e+l 

8 

45 

1.05595 

1.09107e+2 


Table 1: Numerical results for fc-point bilinear quadratures on P n x P n for L 2 on the interior of the 
reference right triangle. 




Figure 1: Evaluation points for bilinear quadratures on P„ x P n for L 2 on the interior of an 
equilateral triangle, for n = 4,8. 


In Figure 1, the evaluation points of two bilinear quadrature rules on the equilateral 
triangle are shown. Notice that the points possess some symmetries. The expectation 
that quadrature points for polynomials should have some symmetries has been exploited 
in the past to reduce the complexity of searching for classical quadratures [XG10]. In the 
quasi-Newton method used to solve (2.8), however, no symmetry conditions were explicitly 
enforced. 


3.4 Numerical accuracy of quadratures on triangles 

In the section the computed bilinear quadratures on triangles are compared against existing 
high-order classical quadrature schemes on triangles in the setting of orthogonal projections. 
Given the space To = P„ on Q with L 2 -orthonormal basis p = {fi}, orthogonal projection 
operator P onto To, and given g G C 00 (n), we wish to compute 


[Pgb 


(/ij3)l 2 


(fk,g)L2 


The column vector [ Pg\p can be computed using either an approximate orthogonal projec¬ 
tion, or a classical quadrature for each entry J fl fig. 

Since the approximate orthogonal projection matrix F*W, weights of the classical quadra¬ 
ture, and locations of evaluation points are all precomputed, the flop cost for each method 
is solely determined by the number of evaluation points needed. The 28-point bilinear 
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P's 

P' 

C 

TP 

Dunavant 

9.38e-14 

3.97e-01 

4.97e-05 

9.06e-03 

Xiao / Gimbutas 

3.29e-15 

2.73e-01 

1.91e-05 

4.74e-03 

Bilinear 

3.92e-15 

3.99e-15 

6.74e-06 

1.71e-03 


Table 2: Average P-norm relative error in computing approximate orthogonal projection coefficients 
onto Pg for four different sets of functions using three methods that use 28 function evaluations. 


quadrature as shown in Figure 1 was utilized. For comparison we chose two different 28- 
point classical quadratures, each exact on polynomials of degree up to 11, due to Dunavant 
[Dun85] and Xiao and Gimbutas [XG10], respectively. These quadratures were computed 
using the libraries available from [Burl5]. Both classical quadratures were similarly trans¬ 
formed to an equilateral triangle of side length 1 . 

For our numerical experiments, we draw the projected function g from four different 
probability distributions of functions, which we denote by Pg,Pg, and TP. 

We define 

K ■= {g e P« ; IIsIIl^ = i}, 

with probability measure given by drawing a random vector of coefficients uniformly in 
[— 1 , l] fc , and then normalizing the coefficients to have £ 2 -norm 1 , and using those as the 
Fourier coefficients on the orthonormal polynomials on the triangle. 

The set C contains smooth functions with slow decay, and is defined by functions of the 
form 

9 ^ v)= l + ( ai x + a 2 y) 2 ’ 

where a = ( 01 , 02 ) is drawn uniformly from the unit circle. 

The set TP contains smooth non-polynomial functions with oscillations, and has ele¬ 
ments of the form 

g{x , y) = e aiX+a2V cos( 4 &! 2 ; + 4 b 2 y)p{x, y), 

where parameters ( 01 , 02 ) and ( 61 , 62 ) are both drawn uniformly from the unit circle, and 
p(x, y) is a random element of P ' 2 with L 2 norm 1 as chosen in the same manner as for the 
first two cases. 

For each randomly chosen function g , we computed the column vector {PQg\p using the 
three quadrature methods. The exact value [Pg\p was computed with a 295-point classical 
quadrature that exactly integrates polynomials up to degree 40, as computed in [XG10]. 
The £ 2 norm relative error was averaged over 10 4 randomly generated g for each of the four 
classes of functions. The resulting average relative errors are shown in Table 2. 

On Pg, all three quadrature rules achieve very high accuracy, with the Dunavant quadra¬ 
ture losing one digit of accuracy and both Xiao/Gimbutas and bilinear quadratures correctly 
computing the orthogonal projection up to double precision. This is expected since all 
quadratures are designed to integrate such polynomial functions exactly. 

On Pg, neither classical quadratures are accurate to full precision because both classical 
quadratures are only capable of exactly integrating polynomials of degree up to 11. Since 
the bilinear quadrature can exactly integrate Pg x P 6 , it has mean error on the order of 
machine precision. 

On the sets C and TP, none of the quadratures are accurate to machine precision since 
none of the functions are polynomials. However, the bilinear quadrature achieves better 
accuracy than the classical quadratures despite having the same number of evaluation points. 

The existing classical quadratures are already very good, integrating non-polynomial 
functions from C and TP with several digits of accuracy. Additionally, the classical quadra¬ 
ture of Xiao/Gimbutas performs better than the Dunavant quadrature in all four cases. 
However, the bilinear quadrature was as good or better than the classical quadratures in 
each case, despite using the same number of evaluations. This result is explained by the fact 
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n 

k 

a 

Koo (F*W) 

0 

1 

0.00000 

1 .00000e+0 

1 

3 

0.67739 

2.91852e+0 

2 

6 

0.79523 

7.50137e+0 

3 

10 

0.92888 

1.17526e+l 

4 

15 

0.97590 

2.68367e+l 

5 

21 

0.99701 

3.14417e+l 

6 

28 

1.00066 

6.42937e+l 

7 

36 

1.00711 

7.34237e+l 

8 

45 

1.00784 

1.03464e+2 


Table 3: Numerical results for fc-point bilinear quadratures on P n x P n for L 2 on the interior of a 
square. 



Figure 2: Evaluation points for bilinear quadratures on P n x P n for L 2 on the interior of a square, 
for n = 4, 8. 


that bilinear quadratures are specifically designed for the orthogonal projection problem, 
while classical quadratures are designed for evaluating a linear functional. 

3.5 Bilinear quadratures on other domains 

In this section bilinear quadratures for L 2 inner products of polynomials on the interiors 
of a square and a circle are computed. We observe that, just as in the case of triangles, 
minimizing according to (2.8) produces well-behaved evaluation points. 

For the case of the square domain [—1, l] 2 , orthogonal polynomials are P n {x)P m {y), 
where P n is the nth Legendre polynomial. Table 3 shows the minimized leading singu¬ 
lar value er and matrix condition number Kqo for several fc-point bilinear quadratures on 
the square. Interestingly, the evaluation points on the square do not appear to obey any 
symmetries. 

Remark. One can produce a classical quadrature scheme on the square by simply taking the 
tensor product of two Gaussian quadratures on an interval. However, this exactly integrates 
basis functions of the form x a yP with 0 < a < n, 0 < f3 < n, rather than integrating 
polynomials whose total degree does not exceed some value. 

On the unit disk, an orthogonal basis of polynomials is given in polar coordinates by the 
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n 

k 

a 

Koo (F*W) 

0 

1 

0.00000 

1 .00000e+0 

1 

3 

0.67617 

3.04857e+0 

2 

6 

0.79868 

5.50559e+0 

3 

10 

0.89712 

1.01509e+l 

4 

15 

0.94133 

1.59179e+l 

5 

21 

0.97804 

2.24193e+l 

6 

28 

1.00337 

3.94055e+l 

7 

36 

1.02908 

5.60579e+l 

8 

45 

1.07413 

6.75064e+l 


Table 4: Numerical results for fc-point bilinear quadratures on P n x P n for L 2 on the unit disk. 




Figure 3: Evaluation points for bilinear quadratures on P n x P„ for L 2 on the unit disk, for n = 4,6. 


Zernike polynomials Z m ^ n {r,d), defined by 

■— Qm,n ) COs(?7i$), Z^ mn ^T,6^ . — Qm,n (f) sin(mfl), 



where n > m > 0 are integers and n — m is even. Table 4 shows the minimized leading 
singular value a and matrix condition number for several fc-point bilinear quadratures 
on the unit disk. 

3.6 Bilinear quadrature for the Sobolev inner product 

In this section we compute bilinear quadratures that evaluate the Sobolev inner product 

(f,g)Hi= f Df(x) ■ A{x)Dg(x) + f{x)g{x)dx, 

Jn 

where A(x) is symmetric positive definite on f 1. One advantage of a bilinear quadrature for 
H 1 is that the above integral can be numerically evaluated using only point evaluations of 
/, g and does not require evaluating any derivatives. 

For = [— 1, 1], bilinear quadratures for H 1 on P ra x P„ and minimal on P n +i l~lP^ were 
computed for two positive weight functions A(x) = l+x 2 and A{x) = e x . Orthogonalization 
was performed by starting with the Legendre polynomials and computing the Gram matrix 
M using a 40-point classical Gaussian quadrature. 
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n 

k 

G 

Koo(F*W) 

i 

2 

0.00000 

5.00000e+0 

2 

3 

0.00000 

1.38132e+l 

3 

4 

0.00000 

2.72011e+l 

4 

5 

0.00000 

6.59254e+l 

5 

6 

0.00000 

1.21461e+2 

6 

7 

0.00000 

1.86818e+2 

7 

8 

0.00000 

2.86549e+2 

8 

9 

0.00000 

4.22824e+2 

Table 5: Numerical results for bilinear quadratures 

on P n xP n for H 

A(x) = 1 + x 2 . 




n 

k 

G 

Koo{F*W) 

i 

2 

0.00000 

4.52560e+0 

2 

3 

0.00000 

1.30446e+l 

3 

4 

0.00000 

2.49183e+l 

4 

5 

0.00000 

5.13338e+l 

5 

6 

0.00000 

9.28987e+l 

6 

7 

0.00000 

1.50063e+2 

7 

8 

0.00000 

2.28284e+2 

8 

9 

0.00000 

3.30651e+2 


1,1] with the weight function 


Table 6: Numerical results for bilinear quadratures on P n xP n for A/ 1 [— 1,1] with the weight function 
A(x) = e x . 


In Tables 5 and 6 the singular value a and condition number Koo are shown for the 
two bilinear quadratures for H 1 . In all cases, cr is zero up to machine precision, since the 
exact solution to the minimization (2.6) is the roots of the (n + l)th-degree ff ^orthogonal 
polynomial, just as for Gaussian quadratures. We observe that the condition number of the 
approximation projection matrix F*W is larger than in the L 2 case. This can be explained 
by the fact that small perturbations in the function values can lead to large perturbations 
in the derivatives. 


4 Conclusions 

A quadrature framework for numerically evaluating a continuous bilinear form on function 
spaces has been presented, and an optimization procedure for computing such quadratures 
has been outlined. We have argued that this is the correct approach to numerically evalu¬ 
ating orthogonal projections of functions onto a fixed subspace. 

We have also observed that the optimization approach for finding bilinear quadratures 
does not depend on the ambient dimension, the domain of integration, or the function space 
to be integrated exactly. Despite this generality, in our numerical experiments we found the 
resulting quadratures perform well, achieving both efficiency and accuracy. 

There are several topics to explore in future work. One is the construction and utiliza¬ 
tion of bilinear quadratures tailored to specific high-order Galerkin methods. Another is the 
investigation of the performance of bilinear quadratures for evaluating other (non-Sobolev) 
bilinear forms. Yet another finding an efficient numerical method for solving the optimiza¬ 
tion problem (2.15) for the non-invertible case. In that case, a bilinear quadrature is not 
uniquely determined by its evaluation points, and the optimization problem gains many 
additional degrees of freedom. Lastly, one could investigate the use of bilinear quadratures 
for solving integral equations. Such quadratures may prove useful in the Nystrom discretiza- 
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tion of Fredholm integral operators [Bol72] or boundary integral equations on domains with 
corners [BRS10]. 
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